#include <stdio.h>
#include <stdlib.h>
#include <math.h>

#define ERR_INF 1e-6	//If |xAccu| is too small, use ERR_INF to alternate

//C Function pointers
typedef double (*FUNC)(double,double *);
typedef double (**pFUNC)(double,double *);
//Function points (t,x1,...,xn)
typedef struct{
	double t;
	double *x;
}POINT,*pPOINT;
//Function solution
typedef struct{
	int n;
	int PointCount;
	pPOINT Curve;	//Solved points as an array
}SOLUTION;
//Solving parameters
typedef struct{
	double a;	//Solving interval: [a,b]
	double b;
	double *x0;	//Initial value at a
	double h;	//Step length
	double eps;	//Relative error
	int n;		//Number of unknown variables
}SOLV_PARAM;

SOLUTION Euler(pFUNC Func,SOLV_PARAM Param);
void ArrayCopy(double *,double *,int);
double VecNorm(double *,int);
void Output(SOLUTION);

double fun1(double t,double *x){
	return x[1];
}
double fun2(double t,double *x){
	return -2*x[0]/(t*t)+2*x[1]/t;
}
int main(){
	double x0[2]={1,4};
	
	SOLV_PARAM Param;
	Param.h = 0.1;
	Param.eps = 1e-6;
	Param.a = 5;
	Param.b = 10;
	Param.x0 = x0;
	Param.n = 2;
	
	pFUNC Func;
	Func = (pFUNC)malloc(Param.n*sizeof(FUNC));
	Func[0] = fun1;
	Func[1] = fun2;
	
	SOLUTION Solution;
	Solution = Euler(Func,Param);
	Output(Solution);
	
	free(Func);
	return 0;
}
//Euler solver core
SOLUTION Euler(pFUNC Func,SOLV_PARAM Param){
	SOLUTION Solution;
	Solution.n = Param.n;
	Solution.PointCount = 0;
	Solution.Curve = NULL;
	
	double t,h,h0;
	double *x,*xMid,*xInAccu,*xAccu,*k1,*k2;
	t = Param.a;
	h = Param.h;
	h0 = Param.h;
	x = (double *)malloc(Param.n*sizeof(double));
	xMid = (double *)malloc(Param.n*sizeof(double));
	xInAccu = (double *)malloc(Param.n*sizeof(double));
	xAccu = (double *)malloc(Param.n*sizeof(double));
	k1 = (double *)malloc(Param.n*sizeof(double));
	k2 = (double *)malloc(Param.n*sizeof(double));
	
	ArrayCopy(x,Param.x0,Param.n);	//Initial value
	
	int i;
	double ErrBase,OutFlag=0;
	while (t<Param.b)
	{
		//The end of the interval
		if ((Param.b-t)<h){
			h = Param.b-t+ERR_INF;
			OutFlag = 1; 
		}
		
		for (i=0;i<Param.n;i++){
			k1[i] = h*(Func[i])(t,x);
			xMid[i] = x[i]+0.5*k1[i];
		}
		//Two integral methods
		for (i=0;i<Param.n;i++){
			k2[i] = h*(Func[i])(t+0.5*h,xMid);
			xInAccu[i] = x[i]+k1[i];				//Euler
			xAccu[i] = x[i]+0.5*k1[i]+0.5*k2[i];	//Euler steplength h/2
		}
		//Calculation error
		for (i=0;i<Param.n;i++)	xMid[i] = xInAccu[i]-xAccu[i];
		
		ErrBase = VecNorm(xAccu,Param.n);
		if (ErrBase<ERR_INF) ErrBase = ERR_INF;
		if (OutFlag == 1)	ErrBase = 5*VecNorm(xMid,Param.n)/Param.eps;
		
		//Step too short
		if ((VecNorm(xMid,Param.n)/ErrBase)<(0.05*Param.eps) && h<h0)
			h *= 2;
		//Step too long
		else if ((VecNorm(xMid,Param.n)/ErrBase) > Param.eps)
			h /= 2;
		//Appropriate step length
		else{
			Solution.PointCount++;
			
			//New point creation
			Solution.Curve = (pPOINT)realloc(
				Solution.Curve,
				Solution.PointCount*sizeof(POINT));
			
			//Assign value to the new point
			t += h;
			(Solution.Curve[Solution.PointCount-1]).t = t;
			(Solution.Curve[Solution.PointCount-1]).x
				= (double *)malloc(Param.n*sizeof(double));
			ArrayCopy((Solution.Curve[Solution.PointCount-1]).x,
				xAccu,Param.n);
			
			//New integral point
			ArrayCopy(x,xAccu,Param.n);
		}
		
		if (OutFlag==1)	break;
		
	}
	
	return Solution;
	
	free(x);
	free(xMid);
	free(xInAccu);
	free(xAccu);
	free(k1);
	free(k2);
}
//Output solution
void Output(SOLUTION Solution){
	int i,j;
	
	FILE *FileOut;
	FileOut = fopen("out.csv","w");
	for (i=0;i<Solution.PointCount;i++){
		fprintf(FileOut,"%f,",Solution.Curve[i].t);
		for (j=0;j<Solution.n;j++){
			fprintf(FileOut,"%f,",(Solution.Curve[i].x)[j]);
		}
		fputc('\n',FileOut);
	}
	
	fclose(FileOut);
}
//Copy array from double* to double*
void ArrayCopy(double *dest,double *src,int n){
	int i;
	for (i=0;i<n;i++)	dest[i] = src[i];
}
//Norm of a double vector
double VecNorm(double *vec,int n){
	int i;
	double Norm=0;
	for (i=0;i<n;i++)	Norm += vec[i]*vec[i];
	return Norm;
}
